﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SHPBaidu
{
    public class BaiduConvert
    {
        private const double PI = 3.14159265358979324;
        private const double EARTH_RADIUS = 6378245.0;
        private const double EE = 0.00669342162296594323;
        private const double X_PI = 3.14159265358979324 * 3000.0 / 180.0;

        /// <summary>
        /// 坐标是否在中国外
        /// </summary>
        /// <param name="lat"></param>
        /// <param name="lng"></param>
        /// <returns></returns>
        private bool OutOfChina(double lat, double lng)
        {
            if (lng < 72.004 || lng > 137.8347) return true;
            if (lat < 0.8293 || lat > 55.8271) return true;
            return false;
        }
        /// <summary>
        /// 纬度转换
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <returns></returns>
        private double TransformLat(double x, double y)
        {
            double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.Sqrt(Math.Abs(x));
            ret += (20.0 * Math.Sin(6.0 * x * PI) + 20.0 * Math.Sin(2.0 * x * PI)) * 2.0 / 3.0;
            ret += (20.0 * Math.Sin(y * PI) + 40.0 * Math.Sin(y / 3.0 * PI)) * 2.0 / 3.0;
            ret += (160.0 * Math.Sin(y / 12.0 * PI) + 320 * Math.Sin(y * PI / 30.0)) * 2.0 / 3.0;
            return ret;
        }
        /// <summary>
        /// 经度转换
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <returns></returns>
        private double TransformLng(double x, double y)
        {
            double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.Sqrt(Math.Abs(x));
            ret += (20.0 * Math.Sin(6.0 * x * PI) + 20.0 * Math.Sin(2.0 * x * PI)) * 2.0 / 3.0;
            ret += (20.0 * Math.Sin(x * PI) + 40.0 * Math.Sin(x / 3.0 * PI)) * 2.0 / 3.0;
            ret += (150.0 * Math.Sin(x / 12.0 * PI) + 300.0 * Math.Sin(x / 30.0 * PI)) * 2.0 / 3.0;
            return ret;
        }
        /// <summary>
        /// 地球坐标转换为火星坐标
        /// </summary>
        /// <param name="wglat">WGS纬度</param>
        /// <param name="wglng">WGS经度</param>
        /// <returns>返回近似火星坐标系</returns>
        public double[] GPS2Mars(double wglat, double wglng)
        {
            if (OutOfChina(wglat, wglng))
                return new double[] { wglat, wglng };
            double dlat = TransformLat(wglng - 105.0, wglat - 35.0);
            double dlng = TransformLng(wglng - 105.0, wglat - 35.0);
            double radlat = wglat / 180.0 * PI;
            double magic = Math.Sin(radlat);
            magic = 1 - EE * magic * magic;
            double sqrtmagic = Math.Sqrt(magic);
            dlat = (dlat * 180.0) / ((EARTH_RADIUS * (1 - EE)) / (magic * sqrtmagic) * PI);
            dlng = (dlng * 180.0) / (EARTH_RADIUS / sqrtmagic * Math.Cos(radlat) * PI);
            double mglat = wglat + dlat;
            double mglng = wglng + dlng;
            return new double[] { mglat, mglng };
        }
        /// <summary>
        /// 采用二分法 火星坐标反算地球坐标
        /// </summary>
        /// <param name="gclat">火星坐标纬度</param>
        /// <param name="gclng">火星坐标经度</param>
        /// <returns></returns>
        public double[] Mars2GPS(double gclat, double gclng)
        {
            double initDelta = 0.01;
            double threshold = 0.000000001;
            double dlat = initDelta;
            double dlng = initDelta;
            double mlat = gclat - dlat;
            double mlng = gclng - dlng;
            double plat = gclat + dlat;
            double plng = gclng + dlng;
            double wgslat = 0;
            double wgslng = 0;
            int i = 0;
            while (true)
            {
                wgslat = (mlat + plat) / 2.0;
                wgslng = (mlng + plng) / 2.0;

                double[] tmps = GPS2Mars(wgslat, wgslng);
                double tmplat = tmps[0];
                double tmplng = tmps[1];

                dlat = tmplat - gclat;
                dlng = tmplng - gclng;

                if ((Math.Abs(dlat) < threshold) && (Math.Abs(dlng) < threshold))
                    break;
                if (dlat > 0) plat = wgslat;
                else mlat = wgslat;
                if (dlng > 0) plng = wgslng;
                else mlng = wgslng;

                i += 1;
                if (i > 10000) break;
            }
            return new double[] { wgslat, wgslng };

        }
        /// <summary>
        /// 火星坐标转换为百度坐标
        /// </summary>
        /// <param name="gglat">火星纬度</param>
        /// <param name="gglng">火星经度</param>
        /// <returns></returns>
        public double[] BD_encrypt(double gglat, double gglng)
        {
            double x = gglng;
            double y = gglat;
            double z = Math.Sqrt(x * x + y * y) + 0.00002 * Math.Sin(y * X_PI);
            double theta = Math.Atan2(y, x) + 0.000003 * Math.Cos(x * X_PI);
            double bdlng = z * Math.Cos(theta) + 0.0065;
            double bdlat = z * Math.Sin(theta) + 0.006;
            return new double[] { bdlat, bdlng };
        }
        /// <summary>
        /// 百度坐标转火星坐标
        /// </summary>
        /// <param name="bdlat"></param>
        /// <param name="bdlng"></param>
        /// <returns></returns>
        public double[] BD_decrypt(double bdlat, double bdlng)
        {
            double x = bdlng - 0.0065;
            double y = bdlat - 0.006;
            double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * X_PI);
            double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * X_PI);
            double gglng = z * Math.Cos(theta);
            double gglat = z * Math.Sin(theta);
            return new double[] { gglat, gglng };
        }

    }
}
